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rapidly  convergent  numerical  method  for  finding  q  is  obtained. 
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SIGNIFICANCE  AND  EXPLANATION 


NUMERICAL  SOLUTION  OF  DEGENERATE  VARIATIONAL  INEQUALITY 
ARISING  IN  THE  FLUID  FLOW  THROUGH  POROUS  MEDIA 


C.  W.  Cryer*  and  s.  z.  Zhou** 

1 •  Introduction 

Moat  steady-state  porous  flow  free  boundary  problems  may  be  reduced  to  elliptic, 
variational  or  quasi-varlational  Inequalities  for  which  the  numerical  solutions  have  been 
studied  by  many  authors  (see,  for  instance,  Baiocchl  and  Capelo  [1978],  Oden  and  KikuChi 
[1979],  and  their  references).  In  some  axisyasntric  problems  there  appears  another  kind  of 
variational  inequalities  -  degenerate  variational  inequalities.  In  this  paper  we  propose  a 
numerical  method  for  a  type  of  degenerate  variational  inequality  arising  in  the 
axisyimaetrlc  porous  flow  well  problems  which  have  been  studied  in  Cryer  and  Zhou  [1961]. 

Ms  use  the  finite  element  method  to  obtain  the  discrete  problem.  The  convergence  of  the 
solution  of  the  discrete  problem  to  the  solution  of  the  degenerate  variational  inequality 
is  proved.  The  solution  of  the  physical  problem  depends  upon  the  unknown  discharge  q.  Ms 
’lve  a  numerical  method  for  finding  q  which  has  been  found  to  oon verge  rapidly. 

Here  ws  recall  some  notations  and  results  about  weighted  Sobolev  spaces  V1  and  v2. 
(See  Chang  and  Jiang  [1978],  Zhou  [I960].) 

A  -  a  bounded  domain  in  (r,z) -plane  with  a  locally  Lips  chits  boundary  T,  and  with 
r  >  0 . 

•  2 

C  (R  >  -  the  space  of  functions  infinitely  differentiable  in  (r,s)-plane. 

C  (A)  «  (v  j  v  is  defined  in  A  and  has  extension  in  C  (R2)) 

Cg(A;T  )  -  (v  6  C  (A)  :  v  «  0  in  some  neighborhood  of  r*),  where  T*  c  p.  if 
T  -  r  then  denote  C0(A|T*)  by  Cg(A) . 
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V^(A)  -  the  closure  of  CQ (A)  in  V^tA). 

1  *  •  *  I 

VQ(A?r  )  -  the  closure  of  CQ(A;r  )  in  V  (A). 

The  following  propositions  can  be  easily  shown. 

Proposition  1.1.  If  T  n  {r  -  o)  -  0  and  meastT  )  >  0  then  there  exists  a  constant  C 
such  that 

Ivl2  <  c  /  ♦  (4~> 2] rdrdx,  V  v  e  v’(A>r  )  . 

V1 (A)  A  3r  »«  0 

Proposition  1 .2 .  v'(A),  V2(A)  and  vJ(A»T  )  are  Banach  spaces,  .here  we  take 

'vl  1  '  ,vl  1  ' 

V0(A,r  ,  V 1  ( A) 

Proposition  1 .3 .  C  (A),  co(A|T  )  are  respectively  dense  in  v'tA),  v’fAiT*), 

Proposition  1.4.  If  T  n  {r  *■  0}  »  0  then  there  exists  a  unique  linear  continuous 
operator  tr  i  v'(A)  *  L2(r  )  such  that  tr  v  •  v  on  T  for  any  v  6  C  (A).  Moreover, 
tr  is  a  compact  operator. 


2 .  Continuous  Problem 

Let  D  be  a  L-shaped  domain  such  as  that  shown  in  Figure  1 ,  Set 
^  »  {(r,z)  :  0  <  r  <  Rg,  0  <  z  <  hg). 

The  original  variational  formulation  of  the  physical  problem  is  as  follows  (Cryer  and 
Zhou  [1901]). 

Problem  (PPW)  .  Let  Rg,  Rj,  hg,  h,,,  H  be  nunbers  such  that  R^  >  Rg  >  0,  H  >  >  hg  >  0 . 

Find  functions  i(r)  and  u(r,z)  such  that 

C°([R0,R11),  w(Rg)  »  hw,  w(",)  -  H  (2.1) 
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r0  -  i(r,z)  :  z  -  RQ  <  r  <  Rt>  , 

ft  *  {(r,z)  :  0  <  z  <  *(r),  RQ  <  r  <  R^)  , 

K  -  (v  e  v1  (tl )  :  v  »  0  on  TUTU  (T  n  3Q>  U  T  U  F  }  . 

I  2  3  4  3  0 


By  using  a  kind  of  Baiocchi  transform 


u  »  u(r,z)  in  n 

■  z  In  D\ff 


w(r,z)  «  /*  (u(r,t>  -  t]dt  in  D  .  (2.8 

Cryer  and  Zhou  [1981]  have  dsrivsd  the  following  degenerate  variational  inequality  with  a 
real  parameter  q. 

Problem  (ppwl ) .  Find  w  6  K  such  that 

-  q  q 

/D  rVwq*V<v  ‘  wq)<Jr<ta 

R  <2.9 

*  %  ’  V  <V  -  "q^Z-h/*  +  V"  -  Vrdr<1*'  W  ’  8  \ 


K  -  (v  e  V1  (D)  s  v  >  0  in  D,  v  <  g  (r,H)  in  D\& .  v  -  g  on  T  }  ,  (2.10) 

q  q  i  q  D 


r°-  “/»  ■ 


gq(r,z)  -  0 


-«e  -r 


H  ,  r 

r +  q  tn  r 


on  r3ur4 


n  ,  »  v 

T  +  q  ln  57  ' 


R„  (h  -z) 


For  Problem  (PPW1 )  there  is  a  equivalent  form  which  is  more  oonvenient  with  regard  to 
the  numerical  solution. 


such  that 


Pro  blent  (PPW2 )  •  Find  w  6  K*‘ 

-  q  q 

J(w  )  -  min  J(v) 

q  .. 

V0K 

q 


where 


J(v)  -  JD  r|Vv|2drd*  -  2(hw  -  hQ)  /Q°  vl^  rdr 


-  2  !  v  rdrdz 
D 


Cryer  and  Zhou  [1961]  have  proved  the  following  results* 


Proposition  2.1.  For  any  q  with  0  <  q  <  q^,  where 


K2  -  (h,  -  hQ)2 

q0  “  2  IntR^Rjj) 


(PPWI )  has  a  unique  solution . 

Proposition  2 .2 .  Let  wq  be  the  solution  of  (PPW1).  Then  there  exist  u.  B  and 

with  a,  B  real  numbers,  6  f  0,  F*  a  continuous  bounded  function,  such  that  the 

following  two  conditions  are  equivalent: 

w  6  C1^})  n  V2  (D) 

q 

F*(q)  -  Bq  -  a  -  0  . 

He  call  Wq  a  regular  solution  of  (PPWI )  if  (2.16)  is  valid. 

Proposition  2.3.  If  Wq  is  a  regular  solution  of  (PPWI)  then  0  <  q  <  q^. 
Proposition  2 .4 .  There  exists  at  least  one  regular  solution  Wq  for  Problem  (ppwi) 

Proposition  2.5.  For  the  regular  solution  wq  of  (PPWI)  we  have 

_  3w_  3w_ 

a>-j3»°,  i„  D 

3w_ 

»  h  -  h.  on  r 
«n  w  0  6 

■  0  on 


-5- 


(2.13) 

(2.14) 


(2.15) 

F*(q) 

(2.16) 

(2. 16’) 


(2.17) 

(2.18) 


Proposition  2.6.  If  u  is  a  solution  of  (PPW)  then  w  defined  by  (2*8)  is  a  regular 


solution  of  ( PPW1 )  corresponding 


rH  9u(r ,t) 

J0  5r 


dt 


Conversely,  if  w_  is  a  regular  solution  of  (PPW1 )  then  the  functions  sP_,  u_  defined  as 

q  q  q 

follows,  is  the  solution  of  (PP'*>: 


^  u  {(r,z)  e  D  :  r  >  rq,  w_  <  g_(r,H)} 

q  q  q 


*  _(r)  -  sup(z  i  (r,z)  e  £!_},  <  r  <  R1 

q  q 


lim  ^_(r) 
r*R,-0  q 


q 


Remark  2.1 .  Physically,  2*q  is  the  "discharge"  while  u  ie  the  "hydraulic  head". 


*  _(R.)  -  lim  *_<r),  f_(R  ) 

q  r"R0f°  q  q 


*  3T1  +  * 

q 


in  D 


3 .  Numerical  Approximation  of  (PPW2 ) 

In  this  section  we  consider  the  approximate  problem  of  Problem  (PPW2>:  Given  q 

and  h  with  0<q<qft,h>0,  find  wh  G  Kh  such  that 
u  q  q 

J(wh)  -  min  J(v)  (3.1) 

q  rflVh 
veK 

q 


where  is  a  convex  closed  nonempty  subset  of  v*1,  and  {v^}  is  a  family  of  finite 


1  h 

dimensional  subspaces  of  V  (D)  .  It  is  not  required  that  K  C  K 

q  q 


So  far,  for  space 


2  j 

V  there  is  no  approximation  theorem  similar  to  that  for  usual  Sobolev  space  H  •  Hence 
we  can  not  prove  the  convergence  theorem  for  our  problem  by  using  usual  methods  such  as 
those  in  Falk  (1974],  Brezi  and  Sacchi  [1976],  Cryer  and  Fetter  [1977]. 
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V 


Let  It  )  be  a  family  of  triangulations  of  the  domain  D.  The  set  of  interior 
n 

gridpoints  will  be  denoted  by  D^,  and  the  set  of  boundary  gridpoints  by  <>D^.  Set 
n  f  .  For  each  Tj,  and  for  each  triangle  T  e  Tj,  we  set 


PCT)  -  diameter  of  T 

0(T)  •  minimise  angle  of  T 

8.  «  minimum  0  (T)  . 

n 


Mb  assume  that 


11m  max  P(T)  *  0  , 

h+0  TOT. 


(3.2) 


and  that  there  exists  a  positive  constant  0Q  independent  of  h  such  that  (Zlamal  [1968]) 


0W  >  an  . 

h  0 


(3.3) 


As  Vh  we  take  the  space  of  linear  finite  elements  corresponding  to  Th  in 
v’jDiI^).  set 

Kh-{w0Vhi»>O  on  D  ,  v  «  g  (r#H)  on  D  -  ft,  ,  v  »  g  on  r  }  .  (3.4) 

q  n  q  n  i  q  on 

He  need  the  following  basic  theorem  (see,  for  instance,  Glowinski  [1980,  Th .  5.2  in 


Ch.  X]  ). 

Ttieoren  3 .1 .  Let  V  be  a  real  Hilbert  space,  a (•  ,•)  -  a  bilinear,  continuous,  symmetric 
and  ooercive  form  on  w  *  v,  f(«)  -  a  continuous,  linear  functional  on  v,  K  -  a  closed, 
convex,  nonempty  subset  of  V,  {X*1}  -  a  family  of  closed,  aonvex,  nonempty  subsets  of  V 
with  Kh  c  vh.  Assume  that 

(i)  If  (vh)  is  bounded  in  V  and  v*1  8  X*1  then  the  weak  cluster  points  of  (vh) 


belong  to  X» 

(ii)  There  exists  a  set  A  C  V 
satisfying  that  vh  8  X*1  and  that 
Then 


with  X 

Hm  vh  • 
h*0 

lim  luh 
h*0 


X  such  that 
strongly  in 

u'v  -  0 


»  v  e  x 
V. 


there  exists 


{vh> 
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where  uh  and  u  are  respectively  the  solutions  of  the  problems 


J(u)  ■  min  J(v) 
v8K 


J(u  )  ”  min  J(v) 
_ h 


provided  that  J(v)  *  a(v,v)  -  2f(v). 

i  —  ** 

In  our  case  we  will  take  C  (D)  n  K  as  the  set  X  in  the  above  theorem. 

q 

From  now  on  we  assume  that  0  <  q  <  q0 . 

lemma  3.2.  C^D)  *  -  K**  in  v’tDI. 

— — —  q  q 

1  —  ** 

Proof.  Set  Y  -  C  (D)  n  Kq  .  By  Cryer  and  Zhou  [1981,  (3.28)]  there  exists  a  function 
vq  e  Y ,  hence  Y  is  nonempty. 

Set 

Y  «  (v  8  v’  (D)  :  v  »  -v  in  0,  v  <  g  <r,H)  -  v  in  d\B  ,  v  -  0  on  T  } 
i  q  q  q  1  D 

y2  -  y,  n  c“<DjrD>  . 

Then  by  the  argument  similar  to  that  of  Lemma  2 .4  in  Glowineki  [I960,  Ch.  II]  we  can  prove 


Y2  ’  *1  (3- 

**  «  * 
v  V  e  we  have  v  •v-vq«Y1.  By  (3.5)  there  exists  a  sequence  vn  e  Y2 

*  *  i  • 

such  that  v„  ♦  v  strongly  in  V  .  Thus  vn  -  vR  +  vq  8  Y  and  vr  *  v  strongly  in 
v1 .  So  we  obtain  that 


1  —  *•  ** 

C  -(D)  OK  D  K 

q  q 


Using  the  well-known  subsequence  argument  and  proposition  1 .4  we  can  easily  obtain 


1  —  »•  •  * 

C  (D)  n  K  D  K 

q  q 


The  following  lemma  may  be  easily  proved. 
Lemma  3.3.  If  T  e  Th  and  vh  8  V*1  then 


/_  v  dr  dz  - 
T 


meas  (T) 
3 


?  h 

l  v («!_>  • 

i-1 


where  mit  (1  -  1,2,3)  are  the  vertices  of  T. 

Now  we  can  prove  the  convergence  theorem* 

Theorem  3 .4 .  Problem  (3.1)  has  a  unique  solution  w*  which  converges  to  the  solution 
wq  of  Problem  (PPW2)  in  V* (D)  as  h  ♦  0. 

Proof .  Let  V  *  V*1,  a(v*,v0)  «  I  r^v  *Vv-  drdz,  and  f(v)  * 

.  *  ~  D  1  2 

R0 

/  rvdrdz  ♦  (h  -  h»)  vl  ^  rdr .  Then  it  is  easy  to  see  that  V  is  a  Hilbert  space 
J  D  w  0  '  0  z«hQ 

with  inner  product 

(Vi,v2)  -  JD  r(v1v2  ♦  VVl-Vv2)drds  , 

that  kh  is  a  closed,  convex,  nonempty  subset  of  V,  that  a(*,M  is  a  symmetric, 
continuous,  coercive  (by  proposition  1*1),  bilinear  form  on  V  x  V,  and  that  f(*)  is  a 
linear,  continuous  (by  proposition  1.4)  functional  on  V.  By  the  well-known  theorem 
(Stampacchia  (1964],  or  Lions  and  Stampacchia  t1967]  we  know  that  Problem  (3.1)  has  a 
unique  solution  wh . 

q 

Now  we  prove  the  convergence  of  w^  by  using  Theorem  3.1.  Let  V  -  V^(D»r^), 

K  -  K**.  It  is  sufficient  to  verify  the  conditions  (1)  and  (il)  because  it  is  obvious  that 

q 

the  rest  of  the  conditions  are  satisfied. 

Verification  of  (1):  Let  fvh}  be  a  sequence  such  that  vh  e  Kh,  and 

-  q 

v11  ♦  v  weakly  in  V^tD*!^)  .  (3.6) 

«* 

We  prove  that  v  e  .  First  we  prove  that 

v  <  g  (r,H)  in  D\S)  .  (3.7) 

q  1 

Define  gq(r,z)  »  9q<Ro**>  in  Let  gh  be  the  plecewisely  linear  interpolation  of 

9q(r,H).  Then  it  is  easy  to  see  that 

gh  ♦  g  (r,H)  in  C°(D)  .  (3.8) 

For  any  ♦  0  Cq(D\Qj)  with  ^  >  0  we  define  ♦  -  0  in  and 

♦h  -  l  t(G  )»(T> 

TOT. 

n 

where  ♦(T)  is  the  characteristic  function  of  T,  and  GT  is  the  centroid  of  T.  It  can 
be  shown  that 

*h  ♦  *  in  C°(D)  .  (3.9) 


*  ‘  )' 
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I 

I 


*  —  _ 

Denote  by  T.  The  union  of  triangles  T  which  are  contained  in  D  «  .  Then  we  have 
n  > 

/«\(j  (vh  *  -  /  .  +  /  . 

1  T.  D\(fl,L'T.  ) 

n  in 


I  L  + 1 


TcT 


DM^uT^ 


I  ♦  (G_)  /,  <vh  "  qNdrd*  +  /  „  (vh  -  gh)*hdrd* 

TcT*  DV(flfTh) 

^  (G  )meas  (T)  3 


I 

TcT 


l  <v (M  )  -  g  (M  ))  +  /  .  (by  lemma  3.3) 

i»1  DVCO^Tj^) 


<  /  t  (v*  -  gh)i(^’drdz  (since  v*1  <  g*1  in  D*1  -  B  ) 

D\(fljjTh) 


(3.10) 


Noting  that  meas(D\(B  U  T*))  *  0  as  h  ♦  0  we  obtain  by  (3.6),  (3.8),  (3.9)  and  (3.10) 
i  n 

that 

/D\n  (V  -  gq(r,H))i|i  drdz  <0,  »  *  0  c'd^B,)  with  *  >  0 

which  implies  that  (3.7)  is  valid. 

By  similar  argument  we  obtain  that 

v  >  0  in  D  .  (3.11) 

Finally,  it  is  easy  to  see  that 

vh  ♦  gq(r,s)  in  cVD)  . 

h  2 

On  the  other  hand,  it  follows  from  proposition  1.4  that  v  ♦  v  strongly  in  L  (f^).  Thus 


we  have 


‘D  * 


(3.12) 


This  equation  as  well  as  (3.11)  and  (3.7)  means  that  v  e  X^. 

1  —  **  _  ** 

Verification  of  (11):  Let  X  »  C  (D)  <">  K  .  By  lemma  3.2  we  have  X  »  K  .  For  any 

1 1  ~ — 1  ■  q  q 

v  0  x  we  take  plecewlsely  linear  interpolation  v*1  of  v.  Then  v*1  0  xjj.  By  a  result  of 
Feng  [1965]  we  have  (note  that  (3.2)  and  (3.3)) 


h  dv  3v  dv  3  v  * 

v  *  v'  TT  *  3r'  5T"  *  57  in  L  (D) 
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V 


q.e.d 


a*  h  ♦  o . 


Therefore  we  have  that 


h 


v 


v  strongly  In  . 


Remark  3.1.  The  condition  (3.3)  may  be  replaced  by  a  weaker  condition  6^  <  *  -  8^,  where 

8*  -  max  8*(T).  8 (T ) •maximum  angle  of  T.  (Cf.  Peng  [1965]  •) 

T«TU 

n 

The  Problem  (3.1)  is  a  quadratic  programming  problem  which  can  be  computed  using 
S.O.R.  with  projection  (Cryer  [1971],  Glowlnski  [1971]).  The  iterative  process  is 
convergent  (see.  for  instance.  Glowinski  et  al.  [1976.  p.  70]).  Me  have  used  Carre's 
scheme  (Carre  [1961])  to  choose  the  relaxation  factor  <•>. 


4.  Numerical  Method  for  Regular  Solution 

Numerical  experiments  indicate  that  w*1  ♦  w  as  h  ♦  0  for  q  with  0  <  q  <  q», 

q  q  ^ 

just  as  theorem  3.4  claims.  But  wq  does  not  satisfy  (2.17)  if  q  does  not  correspond  to 
the  regular  solution.  Here  we  give  a  method  for  searching  for  regular  solution  and 
corresponding  values  of  q. 

The  basic  idea  is  as  follows.  Por  the  solutions  (PPW1 ) .  their  derlvativee  are 
continuous  everywhere  in  D  except  at  the  reentrant  point  P.  The  character  that  a 
regular  solution  possesses  is  that  its  derivatives  are  continuous  even  at  the  point  P. 
Hence  we  may  use  (2.18)  at  P  to  find  regular  solutions.  By  the  way,  equation  (2.40)  in 
Baiocchi  et  al.  [1973]  may  also  be  derived  by  the  same  idea.  Now  we  turn  to  the  concrete 
computation . 

He  choose  a  subset  S  of  the  family  Ct^}  such  that  V  8  S  there  is  a  element 

T*  8  which  has  a  vertical  edge  and  P  is  an  endpoint  of  this  edge.  Denote  by  h 

and  P'  respectively  the  length  and  the  other  endpoint  of  the  edge.  Thus  we  have  the 

discrete  form  of  (2.18)  at  the  point  P: 

f  (q)  -  wh(P)  -  wh(P')  -  h*  (h  -  h.)  -  0  .  (4.1) 

q  q  wo 

This  equation  can  be  numerically  solved  by,  for  example,  the  secant  method  which  is  given 


by 


<*n  '  Vi 


’n+1 


f(v 


"V.1 


f(q_>. 


2,3,, 


(4.2) 
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The  initial  values  q1 ,  q2  must  be  subject  to  the  condition  0  <  q  <  qQ .  Then  we  compute 
Wh  ,  Wh  by  using  S.O.R.  with  projection,  and  ftq,),  f(q2)  by  using  (4.1),  and  q^  by 

q,  q2 

using  (4.2).  This  process  is  repeated  until  prescribed  accuracy  is  reached. 


5.  Numerical  Examples 

We  adopt  the  triangulation  of  D  used  by  Cryer  and  Fetter  [1979]  (see  Figure  2). 
Suppose  that  m  denotes  the  number  of  subdivisions  of  D  in  z-direction,  n  the  number 
of  subdivisions  of  [A2 1  in  r-direction  .  Then  the  coordinates  of  the  gridpoints  are  given 


z -  ( j— 1 )H/m,  1  <  j  <  m  +  1 


r^  -  Rfl  exp[  (i-k-1  )/n*fn(R^/R^ )]  ,  i  »  k+1 ,  . . .  ,n*k+1 


(i-1)R0A,  i  -  1,...,k 


k  -  tR0/Al  , 

4  -  R0t<Ri/Ro>1 


For  vh  e  let  0i:j  ■  ^(r^.z^),  vector  U  *  (u^)  .  Then 


R(i,j) 


JRU,J)(V  ’ 


where  R(i,j)  are  the  rectangles 


R(i,J)  “  (<r,z)  i  ^  <  r  <  ri+1 '  zi  *  *  *  *i+1^  * 


It  is  easy  to  compute 


J„,.  .  v (vh»  -  UTk„(,  .,U  +  2b^  ,  U  . 

R(i, j )  R(i,j)  R(i,j) 


The  matrix  end  the  vector  b^^  j)  are  almost  same  as  that  in  (7.11)  and  (7.12) 

of  Cryer  and  Fetter  [1977].  The  only  difference  is  that  for  j  -  k1  -  1  we  must  add 
respectively  (ho-hw)Ax* xl/2  and  (ho-hw)ix* x^/2 ,  which  correspond  to  the  line  Integral 
in  J(v)  now,  to  b^i  f  j  >  (i ,  J+1 )  and  **r( l ,  j  >  » j+1  > »  where  k1  is  the  value  of  j 

corresponding  z  «  ho . 

For  the  S.O.R.  with  projection  we  only  note  that  there  are  two  constraints  in  our 

problem  -  v*1  *  0  in  (X,  and  v^  4  g  (r,H)  in  D.  -  0 
n  q  n  i 

The  equation  (4.1)  now  becomes 


**£&£*• 


The  results  are  shown  in  the  following  table 


IS- 


Hence 


u(P,)  -  u<P,)  +  0(1-0)P[ 

3  3  dx 


3u(Q0) 

~ 57“' 


By  (A.4)  and  (A. 5)  we  obtain  that 

3v(P)  u<pi>-»<p0>  3U(0) 

"77 - p - 5T"  for  90"*  2  e  P0P, 


3v(P)  u<p2>-“<p3>  PC-P)  .3u(V  *u(V. 

-57 - - tr1  t-H - 57-' 


2s«-L-  <i-> 


*u(Q  )  3u(Q  )  . 

rt9  °oI— 5^ - S7“'  for  *°”®  2  e  P2P3 


Therefore 


.  |i^ai .  i"i£i|  < 


,3vp  .  SulPi,  ,  ,  Wj.  .  MEi,  +  rt,  „o  .  !^ai,  <  „  ♦ 


|v<P)  -  u(P)  I  <  |v(P0)  -  u(P0>|  ♦  I/'  (|I.  *2)dx  ♦  (|l-|H)dy| 


«  (2  +  ctg  a  )»'p  . 


(A. 5) 


“0>“* 


•B  iD  • 
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Appendix  Bi  The  Computer  Program 

implicit  double  precision  (a ,b,c,d,e,f ,h ,1 ,o ,p,q,r ,a ,t ,u,y ) 

integer  ub 

common/ rx/r (150) 

ooomon/sy/s ( 65 ) 

common/ ooefo/cO (150,65) 

oomson/ooef  1/cl  (150,65) 

oommon/ooef  2/c2  (150,65) 

oomaon/coef 3/c3 (150,65) 

oommon/unk/u (150,65) 

common/par/imax,n,m,modit, omega, test, ml ,n1 ,k1 ,k2,kk,q 
.jommn  n/par  1/y  1,  y2  ,y3,b,f  2  ,h2,epa,  iter,  testa 
data  m/8/,n/12/,a/4.8dO/,ab/72.0dO/,y1/48.0dO/, 
q  y2/12.0dO/,y3/9.0dO/,eps/1.0d-6/ 
c  calculate  the  stepsise  for  the  discretisation,  uniformly  for  y, 
c  nonuni  font  ly  for  x 
a1-log(a/a) 
b*a+mb 
b1*log(b/a) 
albl-bl-al 
do  114  ii-1,3 
c  make  the  mesha  finer 
hl-albl/n 
h2-y1/m 

d»a*(exp(h1 )-1 ) 

c  d  is  the  first  dx  in  the  outside  of  the  well 
k-ifix(a/d) 
kk-k+1 
dl-a/k 

c  dl  is  the  dx  for  uniformly  dividing  the  underneath  of 
c  the  well 

n1”n+kk 
ml-m+1 
do  5  i-1,k 

5  r(i)“(i-1)*d1 
do  6  i«kk,n1 

6  r(i)*a»exp( <i-kk)*h1 ) 
do  10  J»1,m1 
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10  S<j)-<j-1)*h2 

c  r(l) ,s( jiare  the  coordinates  of  the  mesh  points 
It  1-dint  (y3/h2)+1 
k2«dint (y2/h2 ) +1 

c  kl  corresponds  the  bottom  of  the  well 
c  k2  correspond  the  water  level 

print  1 S,a ,b,y3 ,y2 ,y 1 ,n ,m,kk ,k1 ,k2 ,eps 
15  format(2x,3hr*^,f12.4,2x,3hre-,f12.4,2x,2hh-,f12.4, 

q  2x, 3hhv-,f 12 .4, 2x, 3hhe-,f 12 .4/ 

q  2x,3hnx-,i5,2x,3hny-,15,2x,3hkk-,i5,2x,3hk1-,15,2x, 

q  3hk2-,15,2x,4heps-,f12.0/) 

c  d  is  now  another  aonstant  (for  saving  storage) 

<Kd1*(y3-y2)/2 

c  calculate  the  term  in  the  right  side  corresponding  the  line 
c  integral 

do  50  1-1 ,k 

c3  ( i,k1 )-  C3(l,k1)+(r(i)+d1/3)»d 
50  c3 ( i+1 ,k1 )-c3 (i+1 ,k1 )+(r (i)+2/3 )+d 

c  calculate  the  coefs  c0,c1,c2  and  the  right  side  term  corres- 
c  ponding  the  multiple  integral 
do  65  1-1 ,n1-1 
do  60  j«1,m 

if  (i  .It.  kk  .and.  j  .gt.  kl-1)  go  to  65 

dx-r(i+1)-r(l) 

dy-s(J+1)-s< j) 

tl-  ( r  ( i )  +dx/3  .Odd )  +dx»dy/2 

t2-dx**2*dy/36.0<J0 

r2-dx*t2+(r ( i) +dx/3 .OdO ) +t1 

r3— *dy*t2/2+(s(  j)  +dy/3 .0d0)*t1 

c0(i, j)-c0(i, j)+t1/dx*»2+t1/dy**2 

c1(i,j)-c1(l,j)-t1/dx«»2 

c2(i,j)-c2(i,j)-t1/dy»+2 

c3(i, j)«c3(i, j  j-t1+(r2-r(i)*t1 )/dx+(r3-s( j )*t1 )dy 

c0 ( i+1 , J )-c0 ( i+1 , j ) +t1/dx**2 

c3(i+1 , j)«c3 (i+1 , j )- (r2-r (i)+t1 )/dx 

c0(  i,J+1  )-c0  (i,  J+1  )+t1/dy**2 

c3(i,j+1)-c3(i,j+1)-(r3-s(j)*t1)/dy 

tip-  (r  ( i )  +2*dx/3  .OdO )  Mx+dy/2 


r2p-dx*t2  +  (r  ( i )  +2*dx/3 .OdO)«t1p 

r3p»- dy»t2/2+(a( j)+2+dy/3 .OdO)+t1p 

cO  (i+1 ,  j+1  )-cO (i+1 , j+1 )+t1p/dx**2+t1p/dy+*2 

c3(i+1,j+1)-c3(i+1,J+1)-t1p-(r2p-r(i+1)*t1p)/dx 

c3 (i+1 , j+1 )-c3 (1+1 , j+1 )-{r3p-a( j+1 )+t1p)/dy 

cO(i,  j+1  )-cO(i,  j+1  )+t1p/dx**2 

c1(t,j+1)-c1(l,j+1)-t1p/dx«*2 

c3 (1, j+1 )-c3(i,  j+1 )+(r2p-r ti+1 )*t1p)/dx 

cO (1+1 , j )-c0 (i+1 , j) +t1p/dy**2 

c2(i+1,J)-c2(i+1,j)-t1p/ay«*2 

60  c3(i+1, j)-c3(i+1,j)+(r3p-a( j+1)«t1p/dy 

65  continue 

c  give  the  firat  initial  value  for  q.  calculate  the  optieuai 
c  acceleration  factor  omega  uaing  carre' a  method 
q-IOO.OdO 
call  lnit 
laax— 150 
oeega-1  .OdO 
call  itera 
omega- 1 .4d0 
do  19  iter-1,20 
call  itera 

if  (iter  .eq.  17)  teat17-teeta 

if  (iter  .eq.  18)  teat18-teeta 

if  (iter  .eq .  19)  teet1>teata 

19  if  (iter  .eq.  20)  teatlO-teata 
pi 8-teat 18/teat 17 
p19*teat19/teet18 
p20-taat20/teat1 9 

if  Up18-p19)*(p19-p20)  .It.  O.OdO)  go  to  20 
if  (dafae(p18-p19)  .1e.  da  be  (pi  9-p20 ) )  go  to  20 
Iamdag-p18-(p1  9-p18)**2/(p18+p20-2*p1  9) 
print  17,laadag 

17  format  (20x,7haitken-,f8.4) 

go  to  25 

20  lamdag-p20 


25 


sq*aqrt{  1  .OdO-  ( laradag+omega- 1  .OdO  )*«2/<  lamdag*oaega**2>) 
omegal-omegaO 
omega0*2  .0d0/(  1  .OdO+sq) 
print  26,omega0 

26  format  (/20x,7homega0», f8 .4 ) 
domega-daba  (omegal  -omegaO ) 

if  ( domega/ ( 2 .OdO -omegaO )  .It.  O.OIdO)  go  to  45 
ome gam -omegaO  -  ( 2 .  OdO  -omegaO  )/4 
print  30,omegam 

30  format  ( 20x ,  Thomegam- ,  f 8 .4 ) 

omega-omegam 
do  40  iter- 1,20 
call  itera 

if  (iter  .eq.  19)  teatOteata 

40  if  (iter  .eq.  20)  teat20-testa 

p20>teat20/teat1 9 
go  to  20 

45  omega-omegaO 
print  46,omega 

46  format  (/2x,6homega-,f8.4) 
teat*1  .OdO 

call  iterat 

f1-f2 

q2*q 

c  give  the  aeaond  initial  value  for  q 
<1*200  .OdO 

c  outer  iteration .  secant  method  for  computing  q 
do  79  lter1»1,10 
call  lnlt 
test»1 .0d0 
call  iterat 
q1-q2 
q2-q 

print  68,lter1,q1,q2,f1,f2 

68  format  (2x,6hiter1-,i3,2x,3hq1-,f14.8,2x,3hq2-,f 14.8,2x, 
q  3hf1-,f14.8,2x,3hf2-,f14.8) 

if  (dabs(f2)  .It.  eps)  go  to  80 
<rq2-(q2-q1)«f2/(f2-f1) 
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c  print  the  final  results 
80  ut*»0 

1b-1 
f-f2 
q-q2 

do  110  n0»1,n1-1.7 
If (nl-1-nO) 115,8S,85 
85  np“min(n1-n0+1 ,7) 

ub-ub+np 

print  90,(r(i) ,i-1b,ub) 

90  format(///1Sx,f1S.8,2x,f15.8,2x,f15.8,2x,f15.8, 

q  2x,f15.8,2x,f15.8,2x,f15.8/) 
do  100  j-1,  a1,1 
jl-nl-j+1 

print  95,  s(jl),  (u(l.jl),  l«1b,ub) 

95  foraat(2x,f 1 1 ,4,2x,f 15.8, 2x,f 15 ,8,2x,f 15,8,2x,f 15.8, 

q  2x,f15.8,2x,f1S.8,2x,f15.8) 

100  oontlnue 

Ib-ub  +  1 
110  oontlnue 

115  print  120,  iterl ,q,f .test 

120  format (///5x,6hiter1«, 13, 2x,2hq-,f 15.8, 

q  2x,2hf-,f15.8,2x,Shtest-,f15.8///> 
do  113  i-t1.nl 
do  113  j-1, ml 
u(i,j)-0.0d0 
c0(i,j)-0 .0d0 
c1(l,j)-0.0d0 
c2(l,j)-0 .0d0 

113  C3(i,  J)-O.OdO 
m-2»m 

114  n-2*n 
stop 
end 

c  Inner  Iteration,  s,o.r.  method  for  computing  u  and  f 
subroutine  Iterat 

Implicit  double  precislon(c,o,t,u,v,q,y,b,f ,h,e) 


oommon/ooefo/cO  ( 97S0 ) 
common/ ooef  1/cl  ( 9750 ) 
common/ cost  2/ c2  (  9750 ) 
ooramon/coef3/c3  (  9750) 
common/unk/u( 9750 ) 

common/par/imax  ,n  ,m  ,modit , omega , test  ,m1 ,01,1(1  ,k2  ,kk,q 

oommon/parl/yl ,y2,y3,b,f2,h2,epa ,itar 

lter-0 

70  iter»iter+1 

moditmaodt  itar ,  10 ) 
if  (modit  .eq.  0)  test-0 .0d0 
do  7  j-2,m 
do  7  i-1 ,n1-1 

if  (i  .la.  kk  .and.  j  .qt.  kl )  go  to  7 

if  (i  .aq.  kk  .and.  j  .aq.  kl)  go  to  7 

ij-i+imax  *  C  j— 1 ) 

im; j-ij-1 

it  -j-ij+1 

ijml-ij-lmax 

ijpl-ij+imax 

c  on  boundary  aagmant  g ama  7,  the  mesh  points  do  not  have 
c  neighbor  mash  points  at  their  laft  aide.  We  use  then  u(1)-0.0d0 
c  Instead  of  u(imlj)  in  the  equations 
if  (i  .eq.  1)  imlj-1 
uold-u(ij) 

unew-  (c3(ij)+c1  (ij)*u(ip1 j )+c2 <ij )*u(ijp1 ) 
q  +c1 (iml j )*u( iml j )+c2 ( ijml )*u( ijml ) )/c0(lj) 
vint“(1 .OdO-omega)*uold-*omega*unet# 
u(i  j  )»dmax1  (vlnt.O  .0d0 ) 

if  (i  .gt.  kk)  u(ij)-*nin1{u(ij),u(i+imax»m)) 
if  (modit  .ne.0)  go  to  7 
vabs- dabs ( u ( 1 J ) -uo Id ) 
testxknaxl  (test.vabe) 

7  continue 

if  (iter  .ge.  500)  go  to  74 
if  (test  .gt.  sps)  go  to  70 

74  print  75, iter, test 

75  format (//20x,5hitsr-,i3,2x,5htest-,f 15.8) 
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i 


i 

n 


c 

c 


f 2-u < kk+i«4K*  (k  1  - 1 ) ) -u ( kk+iiM*«  ( kl -2 )  > -  (y2 _y 3  J >h2 

return 

end 

Initialize  u.  orderj  game  2,  gams  3,  game  4,  gana  5, 
Where  (linear  Interpolation,  and  oonatant  on  gaoa  6) 
subroutine  init 

implicit  double  precision  (b,1,q,r,«,u,y,o,t,f,h) 
ooamon/unk/u (150,65) 
aommon/ rx/r(150) 
ooraaon/ sy/a (65 ) 


elae- 


201 

202 

203 

204 


205 

206 


207 

208 


aontton/par/imax,n  ,m,modit ,  omega ,  teat ,  ml  ,n1  ,k1 ,k2 ,kk,q 

oommon/par  1/y  1  ,y  2  ,y3 ,  b,  f  2  ,  h2 

<to  201  j-1,m1 
u(n1,j)-yl«a(j)-a(j)*»2/2 

*>  202  i-kk.nl 

u(i,a1)-y1»«2/24q«log(r(i)/b) 

*>  203  >k2,m 

u(kk,j)-u(kk,m1) 

*>  204  J-k1,k2 

u(kk,j)«u(kk,m1)-(y2-a(j))**2/2 

do  206  j«2,m 

*>  205  i-kk-f1,nl-i 

lemda-a  ( j  )/y1 

u(i, j)-u(i,a1 )*lamda 

continue 

do  208  j«2,k1 

do  207  i-1,kk 

lamda»B(J)/y3 

u( i « j )“u(kk,k1 )* lamda 

continue 


return 

end 

iteration  for  oarra's  method 
subroutine  itera 
implicit  double  precision 
common/ ooef o/ eo ( 9750) 
common/ roof  1/cl  ( 9750 ) 
common/ ooef  2/c2  (  9750 1 


(«,b,c,d,e,f,h,o,p,q,r ,s,t,u,y) 
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•  y*  :V7^^ppr~ 


non/ooef  3/c3  (  9750 ) 

«on/unk/u  ( 9750 ) 

»on/par/imax,n,ni,modit  .omega,  test  ,m1  ,n1  ,k1  ,k2  ,kk 
■o n/ par 1/y1,y 2, y2, b, f 1 ,h2,eps, iter, testa 
;a-0  .OdO 

ri  j-2,m 
M  i«1,n1-1 

L  .1e.  kk  .and.  j  .gt.  kl)  go  to  71 
L  .eq .  kk  .and.  j  .eq.  kl)  go  to  71 
Umax*( j-1)+i 
|-ij-1 

l«i  j-imax 
l-ij+lmax 
li  .eq.  1)  1*1 j"1 
*-u(ij> 

f-  (c3(ij)+c1(ij)«u(ip1j)+c2<lj)*uUjp1) 

1*1  J  )*u( i*1 j )+c2(ij*1 )*u(ij*1 ) )/c0 (ij ) 
c-(1 .0  dO  -o*ega )  *  no  ld+omega*unew 
))-<knax1  (vint,0  .OdO) 

(1  .gt  •  kk)  u(lj)-daln1 (u(ij) ,u( i+imax**) ) 

*-dabe(u(i  j  )-uold) 

ta-teata+vabs 

>n  that  we  choose  this  norm  for  error  nay  be  found 
[1961] 
irn 
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